Este reporte explora los resultados del monitoreo participativo del estado de salud de árboles de oyamel en el Parque Nacional Desierto de los Leones y sus zonas de influencia en Santa Rosa Xochiac.

Los datos corresponden al resultado del muestreo participativo realizado por 12 brigadistas de Santa Rosa Xochiac, utilizando kobo-conabio como parte del proyecto 308488 Monitoreo y manejo para la conservación de bosques aledaños a la CDMX afectados por contaminación atmosférica de la convocatoria FORDECYT 2019-5.

library(dplyr)
library(tidyr)
library(ggplot2)
library(ggmap)
library(scatterpie)

Los datos del muestreo corresponden a los datos colectados con kobo y limpiados previamente con el script 1_preprocesamiento_datos_kobo.Rmd.

# load data
muestreo_tidy<-read.delim("../data/kobo/muestreo_dic2020_tidy.txt", header = TRUE)
parcelas_tidy<-read.delim("../data/kobo/parcelas_dic2020_tidy.txt", header = TRUE)

# pivot long parcelas data to have health data as a single variable
parcelas_long<-pivot_longer(parcelas_tidy, 
                            cols = healthy:worm, 
                            names_to = "tree_health_simplified",
                            values_to = "n_trees")

Distribución del estado de salud de los árboles por parcela

La siguiente figura muestra el total de árboles muestreados en cada parcela de 10x10 m, y cuántos de estos están bajo alguna categoría de daño:

# Make a nice color pallete and legend order for all plots

my_cols=c("darkgreen", 
              "darkred", 
              "orangered1", 
              "cadetblue", 
              "tan", 
              "beige", 
              "burlywood4", 
              "coral", 
              "aquamarine3", 
              "gray70", 
              "black")

desired_order=c("healthy", 
                "ozone", 
                "ozone_and_other", 
                "others_combined", 
                "drougth", 
                "fungi", 
                "insect", 
                "worm", 
                "acid_rain", 
                "other", 
                "dead")
  
 spanish_labels=c("Sano", 
                  "Ozono", 
                  "Ozono y otros", 
                  "Otros combinados no-ozono", 
                  "Sequía", 
                  "Hongos", 
                  "Insectos", 
                  "Gusano de seda", 
                  "Lluvia acida", 
                  "Otro", 
                  "Muerto")
p <- ggplot(parcelas_long, aes(x=plot, y=n_trees,     fill=tree_health_simplified)) +
  geom_bar(stat="identity") +
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud") 
  

p + theme_bw() +
  ggtitle("Estado de salud de los árboles por parcela") + 
     theme(plot.title = element_text(lineheight=1.1, face="bold")) +
  labs(x="Parcelas", y= "número de árboles") +
  theme(text = element_text(size = 25))

Esta es la distribución de las 48 parcelas:

# code adapted from https://rgraphgallery.blogspot.com/2013/04/rg-plot-pie-over-google-map.html

## configure google api

# You first need to register your api key in https://cloud.google.com/maps-platform/#get-started and follow instructions. The geocoding API is a free service, but you nevertheless need to associate a credit card with the account. Please note that the Google Maps API is not a free service. There is a free allowance of 40,000 calls to the geocoding API per month, and beyond that calls are $0.005 each.
# after you obtain your api, save it in /scripts/api_key.api (not shown in this repo por obvious reasons).

# if you get the following error when running get_map():

#"Error in aperm.default(map, c(2, 1, 3)) : 
#  invalid first argument, must be an array " 

# check this troubleshooting: https://rgraphgallery.blogspot.com/2013/04/rg-plot-pie-over-google-map.html

##  load and register api
api <- readLines("api_key.api")
register_google(key = api)

## plot map
# get map
sat_map = get_map(location = c(lon = -99.3060, lat = 19.2909), zoom = 14, maptype = 'satellite', source = "google")

# plot sampled plots
p_satmap <-  ggmap(sat_map)
p_satmap + geom_point(data=parcelas_tidy,
                      aes(x=X_coordinates_longitude,
                          y=X_coordinates_latitude),
                      color="red") +
          geom_text(data=parcelas_tidy,
                      aes(x=X_coordinates_longitude,
                          y=X_coordinates_latitude,
                          label=plot),
                      color="white",
                 #    check_overlap = TRUE,
                      hjust = 0, vjust=1, nudge_x = 0.0005,
                 size= 5) +
    theme(text = element_text(size = 30))

Según nuestro muestreo, el daño por ozono se distribuye espacialmente de esta forma:

# plot pies in map
p_satmap +
geom_scatterpie(data=parcelas_tidy,
                aes(x=X_coordinates_longitude,
                    y=X_coordinates_latitude,
                    group=plot),
                pie_scale = 1.5,
                cols=desired_order,
                color=NA,
                alpha=1)  +
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud") +
  theme(text = element_text(size = 20))

Haciendo un acercamiento a la zona central:

# zoom plot
sat_map = get_map(location = c(lon = -99.303, lat = 19.2890), zoom = 16, maptype = 'satellite', source = "google")

# plot
p_satmap <-  ggmap(sat_map)

# plot pies in map
p_satmap +
geom_scatterpie(data=parcelas_tidy,
                aes(x=X_coordinates_longitude,
                    y=X_coordinates_latitude,
                    group=plot),
                pie_scale = .7,
                cols=desired_order,
                color=NA,
                alpha=1)  +
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud") +
  theme(text = element_text(size = 30)) +
  
  # add plots id
geom_text(data=parcelas_tidy,
                      aes(x=X_coordinates_longitude,
                          y=X_coordinates_latitude,
                          label=plot),
                      color="white",
                     check_overlap = TRUE,
                      hjust = 0, vjust=1, nudge_x = 0.0005,
                 size= 5)

En total se muestrearon 1765 árboles de los cuales 1722 se encontraban vivos. De estos, el 21.6 % presentó daño por ozono, y el 22.13 presentó daño por ozono en combinación con otro tipo de daño, lo que hace al ozono la fuente de daño más abundante de estos bosques.

A continuación analizamos el daño por ozono en términos de porcentaje de árboles dañados por parcela.

Según la altitud:

# Create new variable with porcentage of ozonoe damage
parcelas_tidy<-parcelas_tidy %>% rowwise() %>% 
                     mutate(., 
                      total=sum(healthy,ozone,ozone_and_other,
                          drougth, acid_rain, other,
                          others_combined, dead, fungi,
                          insect, worm)) %>%
                    mutate(perc.ozone= sum(ozone, ozone_and_other)/total)

#plot
p <- ggplot(parcelas_tidy) +
     geom_point(aes(x=X_coordinates_altitude,
             y=perc.ozone))

p + theme_bw() +
  ggtitle("Daño por ozono según altitud") + 
  labs(x="Altitud de la parcela", 
       y= "Porcentaje de árboles con daño por ozono")+
       theme(plot.title = element_text(lineheight=1.1, face="bold")) +
  theme(text = element_text(size = 16))

O según la latitud. Esto es relevante porque latitudes más al sur están más lejos de la CDMX y por ende de la fuente de contaminantes:

p <- ggplot(parcelas_tidy) +
     geom_point(aes(x=X_coordinates_latitude,
             y=perc.ozone))

p + theme_bw() +
  ggtitle("Daño por ozono según latitud") + 
  labs(x="Latitud de la parcela", 
       y= "Porcentaje de árboles con daño por ozono")+
       theme(plot.title = element_text(lineheight=1.1, face="bold")) +
  theme(text = element_text(size = 16))

Sin embargo, ni la altitud ni la latitud parecen influir en la distribución del daño por ozono.

Distribución del estado de salud de árboles individuales

Examinemos el daño por ozono dependiendo de si la planta fue reforestada o no, y de si se encuentra cubierta o expuesta.

p <- ggplot(muestreo_tidy, aes(x=tree_exposition, 
                      fill=tree_health_simplified)) +
       geom_bar(stat="count") +
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud del árbol") +  
  facet_grid(. ~ reforested, 
             labeller = as_labeller(c("yes" = "reforestado",
                                    "no" = "natural"))) +
  scale_x_discrete(breaks=c("cover", "exposed"),
        labels=c("cubierto", "expuesto"))
p + theme_bw() +
  ggtitle("Estado de salud de los árboles") + 
     theme(plot.title = element_text(lineheight=1.1, face="bold")) +
  labs(x="", y= "Número de árboles") +
  theme(text = element_text(size = 20))

En términos porcentuales:

p <- ggplot(muestreo_tidy, aes(x=tree_exposition, 
                      fill=tree_health_simplified)) +
       geom_bar(stat="count", position = "fill") +
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud del árbol") +  
  facet_grid(. ~ reforested, 
             labeller = as_labeller(c("yes" = "reforestado",
                                    "no" = "natural"))) +
  scale_x_discrete(breaks=c("cover", "exposed"),
        labels=c("cubierto", "expuesto"))
p + theme_bw() +
  ggtitle("Estado de salud de los árboles") + 
     theme(plot.title = element_text(lineheight=1.1, face="bold")) +
  labs(x="", y= "Porcentaje de árboles") +
  theme(text = element_text(size = 20))

Examinemos el daño según la altura de los árboles:

p <- ggplot(muestreo_tidy) +
 scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud del árbol") +
theme_bw()

p2<-p + geom_histogram(aes(x=tree_heigth, 
                      fill=tree_health_simplified))  +
    labs(x="Altura del árbol (m)", y= "Número de árboles") +
    theme(text = element_text(size = 20)) +
    ggtitle("Estado de salud de los árboles según su altura") + 
     theme(plot.title = element_text(lineheight=1.1, face="bold")) 
p2

p2 + 
  facet_grid(. ~ reforested,
    labeller = as_labeller(c("yes" = "reforestado",
                                    "no" = "natural")))

Explorar con mayor detalle los árboles menores a 15 metros:

p <- filter(muestreo_tidy, tree_heigth<15, tree_nodes>0) %>% 
     ggplot(.) +
     scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud del árbol") +
theme_bw()
p + geom_histogram(aes(x=tree_heigth, 
                      fill=tree_health_simplified))  +
    labs(x="Altura del árbol (m)", y= "Número de árboles") +
    theme(text = element_text(size = 20)) +
    ggtitle("Estado de salud de árboles < 15 m según su altura") + 
     theme(plot.title = element_text(lineheight=1.1, face="bold"))

Podemos ver lo mismo, pero con el número de nodos:

pnodos <- p + geom_histogram(aes(x=tree_nodes, 
                      fill=tree_health_simplified))  +
    labs(x="Número de nodos", y= "Número de árboles") +
    theme(text = element_text(size = 20)) +
    ggtitle("Estado de salud de árboles < 15 m \n según su número de nodos") + 
     theme(plot.title = element_text(lineheight=1.1, face="bold"))
pnodos

Y si fueron reforestadas o no:

pnodos + facet_grid(. ~ reforested,
    labeller = as_labeller(c("yes" = "reforestado",
                                    "no" = "natural")))

En términos porcentuales:

pnodos <- p + geom_histogram(aes(x=tree_nodes, 
                      fill=tree_health_simplified),
                      position= "fill", binwidth=1)  +
    labs(x="Número de nodos", y= "Porcentaje de árboles") +
    theme(text = element_text(size = 20)) +
    ggtitle("Estado de salud de árboles < 15 m \n según su número de nodos") + 
     theme(plot.title = element_text(lineheight=1.1, face="bold"))
pnodos

pnodos + facet_grid(. ~ reforested,
    labeller = as_labeller(c("yes" = "reforestado",
                                    "no" = "natural")))
## Warning: Removed 110 rows containing missing values (geom_bar).

También lo podemos explorar por categorías diamétricas:

# changer order of levels of tree_diameter_category
muestreo_tidy$tree_diameter_category<- factor(muestreo_tidy$tree_diameter_category,
levels=c("", "0.5_cm", "1_cm",   "2_cm", "6_cm", "10_cm", "30_cm" , "40_cm")) 

levels(muestreo_tidy$tree_diameter_category)<-c(NA, "0.5_cm", "1_cm",   "2_cm", "6_cm", "10_cm", "30_cm" , "40_cm")

# subset data and plot
filter(muestreo_tidy, tree_heigth<15, !is.na(tree_diameter_category)) %>% 
     ggplot(.) +
     geom_bar(aes(x=tree_diameter_category, 
                      fill=tree_health_simplified),
                      stat="count")  +
    labs(x="Categoría diamétrica", y= "Número de árboles") +
    theme(text = element_text(size = 15)) +
    ggtitle("Estado de salud de los árboles según su diámetro") + 
     theme(plot.title = element_text(lineheight=1.1, face="bold")) + 
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud del árbol") +
    facet_grid(. ~ reforested,
    labeller = as_labeller(c("yes" = "reforestado",
                                    "no" = "natural"))) +
    theme(axis.text.x = element_text(angle=90)) + theme_bw()

En términos porcentules:

# subset data and plot
filter(muestreo_tidy, tree_heigth<15, !is.na(tree_diameter_category)) %>% 
     ggplot(.) +
     geom_bar(aes(x=tree_diameter_category, 
                      fill=tree_health_simplified),
                      stat="count", position = "fill")  +
    labs(x="Categoría diamétrica", y= "Porcentaje de árboles") +
    theme(text = element_text(size = 15)) +
    ggtitle("Estado de salud de los árboles según su diámetro") + 
     theme(plot.title = element_text(lineheight=1.1, face="bold")) + 
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud del árbol") +
    facet_grid(. ~ reforested,
    labeller = as_labeller(c("yes" = "reforestado",
                                    "no" = "natural"))) +
    theme(axis.text.x = element_text(angle=90)) + theme_bw()

Ahora evaluaremos el nivel del daño por ozono, es decir, qué porcentaje del árbol se encuentra dañado, dentro del subset de árboles que están dañados.

my_cols2<-c("gold2", "chocolate1", "orangered", "red4", "darkorchid4")
desired_order_percentage<-c("less than 10%", "10 to 40%", "40 to 50%", "50 to 70%", "more than 70%")

p_od<- muestreo_tidy %>% filter(!is.na(ozone_damage_percentage)) %>%
            ggplot() +
            scale_fill_manual(values= my_cols2, 
                              breaks = desired_order_percentage,
                              labels = c("menos de 10%", "10 a 40%", "40 a 50%",
                                         "50 a 70%", "más de 70%"),
                              name= "Porcentaje del árbol \n dañado por ozono") +
            theme_bw() + theme(text = element_text(size = 20)) 

p_od +
  geom_bar(aes(x=tree_diameter_category,
               fill=ozone_damage_percentage)) +
  labs(x="Categoría diamétrica", y= "Número de árboles") +
  ggtitle("Porcentaje del árbol dañado en árboles con daño por ozono") + 
  theme(plot.title = element_text(lineheight=1.1, face="bold"))

Viéndolo por la edad de los árboles (en árboles <15 años):

p_od +
  geom_bar(aes(x=tree_nodes,
               fill=ozone_damage_percentage)) +
  labs(x="Número de nodos", y= "Número de árboles") +
  ggtitle("Porcentaje del árbol dañado por ozono en árboles <15 m") + 
  theme(plot.title = element_text(lineheight=1.1, face="bold"))
## Warning: Removed 102 rows containing non-finite values (stat_count).

En términos porcentuales:

p_od +
  geom_bar(aes(x=tree_nodes,
               fill=ozone_damage_percentage),
           position = "fill") +
  labs(x="Número de nodos", y= "Porcentaje de árboles") +
  ggtitle("Porcentaje del árbol dañado por ozono en árboles <15 m") + 
  theme(plot.title = element_text(lineheight=1.1, face="bold"))
## Warning: Removed 102 rows containing non-finite values (stat_count).

Y dividido por reforestadas y naturales:

p_od +
  geom_bar(aes(x=tree_nodes,
               fill=ozone_damage_percentage),
           position="fill") +
  labs(x="Número de nodos", y= "Número de árboles") +
  ggtitle("Porcentaje del árbol dañado por ozono en árboles <15 m") + 
  theme(plot.title = element_text(lineheight=1.1, face="bold")) +
  facet_grid(. ~ reforested, 
             labeller = as_labeller(c("yes" = "reforestado",
                                    "no" = "natural")))
## Warning: Removed 102 rows containing non-finite values (stat_count).

p_od +
  geom_bar(aes(x=tree_nodes,
               fill=ozone_damage_percentage),
           position="fill") +
  labs(x="Número de nodos", y= "Porcentaje de árboles") +
  ggtitle("Porcentaje del árbol dañado por ozono en árboles <15 m") + 
  theme(plot.title = element_text(lineheight=1.1, face="bold")) +
  facet_grid(. ~ tree_health_simplified, 
             labeller = as_labeller(c("ozone" = "Sólo daño por ozono",
                                    "ozone_and_other" = "Daño por ozono y otro")))
## Warning: Removed 102 rows containing non-finite values (stat_count).